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Several  papers  on  ionic  liquids  have  been  published  or  submitted  as  a  result  of  this  grant.  The  first  set  of  papers  used  second  order 
perturbation  theory  to  study  the  geometries  and  substituent  effects  of  the  cations  commonly  employed  in  energetic  ionic  liquids.  These  are 
variously  substituted  tria2olium,  tertazolium,  and  pentazolium  cations.  The  heats  of  formation  of  all  species  were  predicted  using  G2  and 
G3  theory.  It  was  consistently  found  that  the  most  energetic  substituents  are  -CN  and  -N3.  When  one  cation  is  combined  with  one  anion, 
proton  transfer  almost  always  occurs  with  no  intervening  energy  barrier,  yielding  a  neutral  pair.  When  two  ion  pairs  are  considered, 
conceptually  similar  to  the  face  of  a  crystal,  the  ion  separated  species  is  predicted  to  be  6  kcal/mol  lower  in  energy  than  the  double  proton 
transferred  neutral  species.  These  calculations  were  done  using  coupled  cluster  theory,  but  this  level  of  theory  is  too  expensive  to  study 
larger  clusters.  Therefore,  we  have  turned  to  the  fragment  molecular  orbital  method,  which  can  retain  the  accuracy  of  full  electronic 
structure  theory  calculations,  while  greatly  reducing  the  cost. 
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This  is  the  final  report  for  the  High  Energy  Density  Matter  (HEDM)  grant,  FA9550-06-1- 
0043,  “Potential  Energy  Surfaces  and  Dynamics  of  High  Energy  Species’"  that  ended  in 
Fall  2008. 

Ionic  Liquids.  Several  papers  on  ionic  liquids  have  been  published  or  are  in  press  as  a 
result  of  this  grant: 

E  Michael  W.  Schmidt,  Mark  S.  Gordon,  and  Jerry  A.  Boatz,  “Triazolium- Based 
Energetic  Ionic  Liquids”,  J.  Phys.  Chem.,  A 109,  7285  (2005). 

The  results  presented  in  this  paper  address  several  different  aspects  of  energetic  ionic 
liquids.  The  conclusions  reached  include  the  following: 

A.  The  contribution  of  the  two  resonance  structures  in  the  unsubstituted  triazolium  cation 
is  roughly  equal,  with  approximately  3/2  p  electrons  found  on  both  charged  nitrogens. 
Substitution  at  ring  nitrogen  atoms  produces  little  change  in  the  extent  of  this  charge 
delocalization,  as  shown  by  both  very  small  geometry  changes  in  the  rings,  and  by 
MCSCF  analysis  of  the  resonance  hybrids.  The  two  substituents  that  cause  the  most 
degradation  to  resonance  are  CN  and  NF2,  but  to  only  a  small  extent.  It  appears  that 
substitution  will  not  hinder  the  ability  of  these  cations  to  form  ionic  liquids. 

B.  A  number  of  energetic  and  less  energetic  substituents  at  the  ring  nitrogens  were 
investigated,  and  in  all  cases  1,4  substituted  triazoliums  are  lower  in  energy  than  1 ,2 
substituted  triazoliums.  There  is  a  small  energy  preference  for  monosubstilution  at  N] 
over  N4  in  the  1 ,4  substituted  triazoliums.  The  heats  of  formation  for  substituted  gas 
phase  triazolium  cations  were  presented,  and  the  CN  substituent  is  suggested  as  a 
synthetic  target,  over  N3,  since  the  former  does  not  require  creating  a  ring  exo  NN  bond, 
but  rather  a  NC  bond. 

C.  Triazolium  dinitramide  systems  are  rich  in  structure.  The  numerous  geometries  and 
their  relative  energies  are  understandable  by  consideration  of  hydrogen  bonds  and 
protonation  energetics.  The  presence  of  only  small  energy  barriers,  or  often,  spontaneous 
transfer  of  protons  from  ionic  dimers  to  produce  neutral  pairs  implies  that  depro  to  nation 
is  a  fundamental  mechanism  for  triazolium  decomposition. 

2.  Deborah  D.  Zorn,  Jerry  A.  Boatz  and  Mark  S.  Gordon,  “Tetrazolium- Based  Energetic 
Ionic  Liquids”,  J.  Phys.  Chem.,  BEK).  1 1 1 10  (2006). 

This  work  followed  the  general  approach  of  that  in  the  previous  paper.  The  cyclic 
structures  are  generally  the  lowest  in  energy  for  this  species;  however,  the  acyclic  for  is 
generally  not  very  much  higher  in  energy.  The  DFT  and  MP2  energies  are  generally  in 
good  agreement  with  each  other.  An  MCSCF  it  orbital  analysis  indicates  that  the 
electrons  in  the  cation  ring  are  delocalized.  Calculated  heats  of  formation  show  that  the 
tetrazolium  cation  ring  has  the  potential  to  release  large  amounts  of  energy  during 
decomposition  and  thus  has  excellent  potential  as  a  high  energy  fuel.  This  is  especially 
true  when  the  ring  is  substituted  with  -N3  or  -CN.  When  a  cation  is  paired  with  oxygen 
rich  anions,  a  single  gas  phase  ion  pair  was  not  generally  found  to  be  stable.  A  proton 
transfers  without  barrier  from  the  cation  to  the  anion  to  form  a  neutral  pair. 


3.  Ian  S.O  Pimienta,  Sherrie  Elzey,  Jerry  A.  Boatz  and  Mark  S.  Gordon,  “Pentazole- 
Based  Energetic  Ionic  Liquids:  A  Computational  Study”,  J.  Phys.  Chem.  A,  J_U,  691 
(2007). 

The  pentazole  cation  has  nitrogen  at  all  positions  on  the  five-member  ring.  Ionic  liquid 
dimers  composed  of  the  five-membered  monosubstituted  pentazole  cation  and  the 
dinitramide,  nitrate,  and  perchlorate  anions  were  investigated  in  this  work.  Of  the  two 
monosubstituted  pentazole  cation  isomers,  it  is  predicted  that  the  1,3-isomer  is  lower  in 
energy  than  the  1,2-isomer.  The  relative  activation  energies  and  decomposition  energies 
for  the  decomposition  of  these  cations  to  N2  +  the  corresponding  azidinium  cation  can  be 
related  to  the  electron-donating  character  of  the  substituent.  An  electron-donating 
substituent  tends  to  increase  the  activation  energy  and  decrease  the  decomposition 
energy.  The  MCSCF  71  orbital  analysis  suggests  that  the  electrons  in  the  cationic 
pentazole  ring  are  delocalized.  The  calculated  enthalpies  of  formation  show  that  electron- 
withdrawing  groups  such  as  -CN  and  -N3  help  to  increase  the  heats  of  reaction.  Several 
dimer  structures  composed  of  the  H-substituted  pentazole  and  either  of  the  anions 
mentioned  above  are  proposed  and  found  to  be  very  unstable.  The  low  barrier  connecting 
the  ionic  and  neutral  dimers,  and  the  large  corresponding  exothermicity  suggests  that  a 
proton  is  spontaneously  transferred  to  the  anion.  So,  a  fundamental  step  for  the 
decomposition  of  protonated  pentazole  is  likely  to  be  deprotonation.  However,  the 
introduction  of  additional  ion  pairs  are  likely  to  stabilize  the  charge  separation  in  the 
ionic  species.  It  is  predicted  that  the  1  -H,2-H-pentazolium-nitrate  pair  with  a  barrier  for 
proton  transfer  of  4.7  keal/mo!  is  the  most  stable  dimer  amongst  all  of  the  ion-pairs 
studied.  The  ion-pair  comprised  of  the  pentazolium  cation  and  perchlorate  anion 
possesses  the  largest  energy  content  and  would  have  been  a  desirable  high-energy  ionic 
liquid  but  the  2.7  keal/mol  barrier  for  proton  transfer  is  prohibitively  small  for  the 
synthesis  and  detection  of  this  ionic  liquid.  The  l-H,3-H-pentazolium-dinitramide  pair  is 
more  ideal  with  a  relatively  larger  barrier  of  4.5  keal/mol  and  a  large  energy  content. 
Among  all  the  ion-pairs  studied,  it  is  calculated  that  the  pentazolium-nitrate  pairs  provide 
the  largest  energy  released  with  a  heat  of  reaction  of  about  130  keal/mol. 

4.  Hui  Li,  Jerry  A.  Boatz,  and  Mark  S.  Gordon,  “Cation-cation  TI-TI  stacking  in  small 
ionic  clusters  of  1,2,4-triazolium”,  J.  Am.  Chem.  Soc.,  130,  392  (2008). 

In  this  paper,  a  novel  combined  MP2/CCSD(T)  method  was  used  to  illustrate  that 
although,  as  noted  above,  single  cation/anion  pairs  tend  to  spontaneously  convert  to  the 
corresponding  neutral  pairs,  double  cation/anion  pairs  are  significantly  more  stable  than 
the  corresponding  neutral  pairs.  In  addition,  ab  initio  calculations  suggest  that  cation- 
cation  7r-7r  stacking  structures  can  exist  in  very  small  ionic  clusters  such  as  two  1 ,2,4- 
triazolium  cations  and  two  dinitramide  or  chloride  anions.  The  structure  motifs  and 
interaction  patterns  provide  new  understanding  of  ionic  materials  with  aromatic  rings.  In 
part  icular,  the  presence  of  two  ion  pairs  appears  to  function  like  a  face  of  a  crystal, 
resulting  in  additional  stability. 


5.  Mark  S.  Gordon,  Jonathan  M.  Mullin,  Spencer  R.  Pruitt,  Luke  B.  Roskop,  Lyudmila  V. 
Slipchenko  and  JerryA.  Boatz,  “Accurate  Methods  for  Large  Molecular  Systems”,  J. 

Phys.  Chem.  (Invited  Centennial  Feature  Article),  in  press 

Three  exciting  new  methods  that  address  the  accurate  prediction  of  processes  and 
properties  of  large  molecular  systems  are  discussed.  The  systematic  fragmentaton  method 
(SFM)  and  the  fragment  molecular  orbital  (FMO)  method  both  decompose  a  iarge 
molecular  system  (e.g.,  protein,  liquid,  zeolite)  into  small  subunits  (fragments)  in  very 
different  ways  that  are  both  designed  to  retain  the  high  accuracy  of  the  chosen  quantum 
mechanical  level  of  theory  while  greatly  reducing  the  demands  on  computational  time 
and  resources.  Each  of  these  methods  is  inherently  scalable  and  is  therefore  eminently 
capable  of  taking  advantage  of  massively  parallel  computer  hardware,  while  retaining  the 
accuracy  of  the  corresponding  electronic  structure  method  from  which  it  is  derived.  The 
effective  fragment  potential  (EFP)  method  is  a  sophisticated  approach  for  the  prediction 
of  non-bonded  and  intcrmolecular  interactions.  Therefore,  the  EFP  method  provides  a 
way  to  further  reduce  the  computational  effort  while  retaining  accuracy,  by  treating  the 
far  field  interactions  in  place  of  the  full  electronic  structure  method.  The  performance  of 
the  methods  is  demonstrated  using  applications  to  several  systems,  including  benzene 
dimer,  small  organic  species,  pieces  of  the  alpha  helix,  water,  and  ionic  liquids. 

Cryogenic  Species.  Timothy  J.  Dudley  and  Mark  S.  Gordon,  “Theoretical  Study  of  the 
Formation  and  Isomerization  of  ALFL”,  Mol.  Phys.,  104,  751  (2006). 

There  has  previously  been  considerable  interest  in  embedding  small,  light  metals  (such  as 
Al)  in  solid  hydrogen.  In  this  work,  the  Iowest-energy  singlet  and  triplet  states  of  the  title 
molecules  were  studied  extensively.  This  is  the  first  time  that  all  of  the  minima  on  the 
two  surfaces  have  been  characterized  at  the  same,  high  level  of  theory  -  multi-reference 
perturbation  theory.  Although  the  singlet  surface  is  generally  lower  than  the  triplet 
surface  at  the  detected  stationary  point,  the  energy  difference  between  the  two  states  is 
less  than  20  kcal/mol.  The  global  minimum  on  the  singlet  potential  energy  surface  is  the 
di-bridged  isomer,  with  other  singlet  isomers  lying  slightly  higher  in  energy.  A  purely 
attractive  singlet  reaction  channel  involving  the  insertion  of  H2  into  the  At-Al  bond  of  an 
excited  Al2  species  is  predicted  to  be  exothermic  by  ~40  keal/moi. 

MCSCF  Hessians.  Timothy  J.  Dudley,  Ryan  M.  Olson,  Michael  W.  Schmidt,  and  Mark 
S.  Gordon,  “Parallel  Coupled  Perturbed  CASSCF  Equations  and  Analytic  CASSCF 
Second  Derivatives”,  J.  Comp.  Chem.,  27,  353  (2006). 

Because  MCSCF  wave  functions  are  very  important  for  many  applications,  as  described 
in  the  previous  paragraph.  In  order  to  characterize  electronic  states  at  this  level  of  theory, 
it  is  important  to  have  analytic  second  derivatives,  Hessians.  MCSCF  Hessians  are  also  a 
key  first  step  in  the  implementation  of  non-adiabatic  coupling  (vibronic  coupling)  matrix 
elements  which  are  very  important  for  the  study  of  surface  crossings  (e.g.,  in 
photochemistry).  Since  MCSCF  calculations  can  be  very  resource  (e.g.,  time,  memory, 
disk)  consuming,  it  is  important  to  try  to  reduce  these  costs.  This  was  accomplished  by 


developing  and  implementing  the  MCSCF  Hessians  as  a  parallel  code.  Good  scalability 
up  to  128  processors  was  demonstrated. 

GAMESS.  Mark  S.  Gordon  and  Michael  W.  Schmidt,  “Advances  in  Electronic  Structure 
Theory:  GAMESS  a  Decade  Later”,  Theory  and  Applications  of  Computational 
Chemistry,  Ch..  41,  C.  E.  Dykstra,  G.  Frenking,  K.S.  Kim,  G.E.  Scuseria,  Eds.,  Elsevier, 
2005. 

The  electronic  structure  code  GAMESS  (General  Atomic  and  Molecular  Electronic 
Structure  System)  has  been  supported  by  AFOSR  for  ~25  years.  In  this  paper,  a 
description  of  new  developments  in  GAMESS  for  the  preceding  10  years  is  presented. 
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PROJECT  SUMMARY/ABSTRACT 
POTENTIAL  ENERGY  SURFACES  AND  DYNAMICS  OF 
HIGH  ENERGY  SPECIES 


The  central  goals  of  the  proposed  research  are  to  (a)  develop  efficient  and  effective 
methods  for  electronic  structure  theory  calculations  and  for  interfacing  such  calculations 
with  dynamics  and  condensed  phase  models,  and  (b)  apply  these  methods  to  design  and 
evaluate  potentially  new  high  energy  species.  The  thermodynamic  properties  to  be 
predicted  are  relative  isomer  energies,  heats  of  formation,  and  ultimately  the  specific 
impulse  (Lp).  The  dynamical  properties  to  be  investigated  include  transition  state 
structures  and  vibrational  frequencies,  classical  barrier  heights,  activation  energies  and 
free  energies,  minimum  energy  reaction  paths  (MEP's),  vibrational  frequencies  and 
curvature  along  the  MEP,  coupling  of  the  transverse  vibrations  with  the  reaction  path  and 
with  each  other,  as  well  as  coupling  with  other  electronic  potential  energy  surfaces.  The 
effect  of  solvation  (e.g.,  exposure  to  water)  and  other  environmental  impacts  will  be 
studied  with  the  aid  of  the  effective  fragment  potential  (EFP)  model.  The  specific 
systems  to  be  explored  include  ionic  liquids  and  a  range  of  high-nitrogen  content  and 
nitrogen-oxygen  content  species.  Polyhedral  oligomeric  silisesquioxanes  are  also  of 
interest  as  protective  coatings  that  are  resistant  to  extreme  environments.  Many  of  the 
proposed  studies  will  involve  collaborations  with  colleagues  at  the  Air  Force  Research 
Laboratory  and  with  other  HEDM  researchers. 

An  essential  ingredient  in  achieving  these  goals  is  our  continuing  development  of 
parallel  electronic  structure  methods.  As  implemented  in  the  electronic  structure  code 
GAMESS,  the  speedups  attained  to  date  are  impressive.  Many  of  the  results  reported  in 
the  previous  results  section  would  have  essentially  been  impossible  without  the 
availability  of  parallel  GAMESS  and  the  HPC  Centers  operated  by  the  HPCMO.  Further 
development  of  parallel  methods  for  correlated  wave  functions  (e.g.,  energy  derivatives 
for  general  and  state-averaged  multi-configurational  (MCSCF)  and  configuration 
interaction  (Cl)  wave  functions,  open  shell  second  order  perturbation  theory),  as  well  as 
related  derivative  coupling  matrix  elements  will  be  very  important.  Likewise,  energy 
derivatives  for  the  generalized  effective  fragment  potential  method  and  interfaces 
between  EFP  and  the  high-level  electronic  structure  methods  will  enable  the  study  of 
environmental  effects  on  HEDM  species  in  both  ground  and  excited  electronic  states. 
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STATEMENT  OF  OBJECTIVES 


We  will  initiate  a  broad  new  effort  in  the  field  of  energetic  ionic  liquids  (EIL).  We  will  focus 
on  both  sophisticated  electronic  structure  calculations  and  on  the  use  of  our  effective  fragment 
potential  (EFP)  model  in  order  to  bridge  the  gap  between  the  gas  and  condensed  phase. 

Properties  of  cation  monomers.  The  first  step  will  be  to  systematically  predict  the  G3  proton 
affinities  for  simple  amines,  R-NH2  +  H*  NH3R+.  This  simple  reaction  will  provide  an  initial 
indication  about  the  inherent  ability  of  a  particular  substituent  to  increase  the  stability  of  the 
amino  proton  and  therefore  resist  proton  transfer  to  the  companion  anion.  The  next  step  will  be 
to  optimize  the  MP2  geometries  of  the  corresponding  substituted  triazoles  and  tetrazoles, 
followed  by  the  prediction  of  the  G3  heats  of  formation.  These  calculations  will  help  identify 
those  substituents  that  are  most  likely  to  increase  the  heat  of  formation  relative  to  the  parent 
compound.  For  those  substituents  that  are  the  most  promising,  the  corresponding  di¬ 
substitutions  will  be  investigated  as  well.  Although  the  main  interest  is  in  substitution  at  N,  the 
same  substituents  will  be  considered  at  C  as  well. 

Cation-anion  pair  interactions.  Once  the  inherent  heats  of  formation  as  a  function  of 
substitution  are  understood,  the  next  step  will  be  to  investigate  the  ion  pair  monomers. 
Therefore,  the  single  ion  pair  calculations  described  above  will  be  followed  by  dimer 
calculations  to  determine  whether  the  presence  of  two  cation-anion  pairs  serves  to  stabilize  the 
species  by,  for  example,  inhibiting  proton  transfer  or  charge  transfer.  Analogous  calculations  to 
those  described  above  for  the  single  ion  pairs  will  be  repeated  for  the  dimer  pairs.  Ultimately, 
the  goal  is  to  build  a  unit  cell  for  each  of  the  most  promising  species,  based  on  their  calculated 
thermodynamic  properties.  The  potential  energy  surfaces  of  these  species  will  provide  a  first 
(ab  initio )  step  in  the  study  of  the  condensed  phases  of  these  systems 

Decomposition  and  oxidation  mechanisms.  All  of  the  foregoing  analyses  focus  on  structural 
predictions  and  thermodynamic  quantities  for  the  parent  compounds  and  ion  pairs.  The  next 
steps  will  be  to  initiate  investigations  of  the  decomposition  and  oxidation  mechanisms  of  the 
cationic  and  neutral  triazole  and  tetrazole  rings. 

Extended  systems  and  connection  with  condensed  phase.  We  will  use  very  high  levels  of  theory 
to  develop  accurate  potential  energy  surfaces  that  can  then  be  used  by  condensed  phase  theorists 
to  build  reliable  potentials  with  which  to  perform  their  simulations.  The  proposed  calculations 
that  are  described  above  are  expected  to  serve  in  this  capacity,  in  addition  to  providing 
important  insights  in  their  own  right.  Second,  in  order  to  build  potentials  that  are  as  reliable  as 
possible,  we  will  develop  new  efficient,  highly  scalable  methods,  so  that  the  most  sophisticated 
methods  are  applicable  to  realistic  problems  on  tractable  time  scales.  Third,  this  group  has,  over 
the  past  several  years,  been  developing  a  very  sophisticated  model  potential,  the  effective 
fragment  potential  (EFP)  that  is  firmly  based  in  quantum  mechanics.  The  EFP  method,  when 
completed,  will  contain  all  of  the  fundamental  physics  involved  in  intermolecular,  and  inter¬ 
ionic,  interactions.  Combined  with  periodic  boundary  conditions  and  either  FMM  or  Ewald 
sums,  the  EFP  method  will  allow  us  to  simulate  the  crystal  phase  changes,  the  melting  process 
and  subsequent  condensed  phase  behavior,  using  a  combination  of  Monte  Carlo  and  molecular 
dynamics  methods. 
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I.  Results  from  Previous  Grant  Period 


A.  THEORY/MODEL  DEVELOPMENT 

Parallel  Developments.  A  highly  scalable  code  for  frozen  core  second  order  perturbation 
theory  gradients  for  closed  shell  molecules  has  been  developed'  and  is  on-line  for  general 
use  at  several  IBM  SP2  and  Cray  T3E  systems,  including  the  T3Es  at  ERDC,  AHPCRC 
and  NAVO  and  the  SP2  at  MHPCC.  This  new  code  uses  the  Distributed  Data  Interface 
(DDI)2,  so  that  the  large  arrays  do  not  have  to  be  replicated.  This  means  that  much  larger 
molecules  can  be  studied  using  geometries  based  on  correlated  wavefunctions.  The  code 
scales  very  well  for  as  many  as  5 12  nodes  and  has  already  been  applied  to  several 
challenging  compounds.  The  development  of  analogous  codes  for  molecules  with 
unpaired  electrons,  using  both  spin-restricted  and  unrestricted  wavefunctions,  is  in 
progress.  The  unrestricted  second  order  perturbation  theory  parallel  code  under  DDI  has 
been  completed3,4  and  is  now  on-line  at  the  same  centers.  A  paper  on  the  restricted  open- 
shell  derivation  has  been  published5,  and  the  development  of  the  parallel  code  is  in 
progress.  A  paper  on  a  DDI  implementation  of  the  SCF  method  has  also  been  published6. 
Improved  parallel  MCSCF  codes  are  under  development,  and  a  distributed  parallel  full  Cl 
code  has  been  implemented  and  a  paper  published7.  The  multi-reference  perturbation 
theory  code  (MRMP2)  has  now  been  implemented  under  DDI.  The  effective  fragment 
potential  (EFP)  code  has  been  made  scalable  under  DDI8.  In  an  important  related 
development,  parallel  GAMESS  now  runs  on  clusters  of  PCs  and  Macs  (running  Linux) 
and  high-end  workstations.  This  is  more  challenging  than  self-contained  massively 
parallel  computers,  since  the  overhead  due  to  inter-node  communications  is  more 
complex.  This  issue  has  been  solved  to  some  extent  by  using  a  Gigabit  Ethernet  switch 
with  large  data  packets.  We  are  also  exploring  alternative  communications  solutions, 
such  as  Myrinet,  SCI,  and  Infiniband.  These  developments  and  their  applications  have 
been  enhanced  by  the  construction  of  a  new  32-node,  128-CPU  IBM  Power3+  cluster 
using  funds  provided  by  a  DURIP  grant.  As  for  all  other  GAMESS  developments,  we 
will  make  our  experiences  in  developing  scalable  clusters  available  to  all  users. 

Condensed  Phase  Methods.  We  have  already  shown  that  our  effective  fragment  potential 
(EFP)  method  for  solvation  is  excellent  for  water,  in  a  variety  of  applications9.  We  are 
now  working  on  extending  the  capabilities  of  the  method  in  several  ways.  We  are 
exploring  several  alternative  approaches  for  incorporating  dispersion  and  other  higher 
order  terms  into  the  method.  Such  terms  are  particularly  important  for  nonpolar  solvents. 
We  are  also  in  the  process  of  extending  the  model  so  that  it  is  equally  applicable  and 
accurate  for  any  solvent.  Key  to  the  success  will  be  the  derivation  of  general  expressions 
for  charge  transfer  and  dispersion  contributions  that  contain  no  fitted  parameters.  This 
will  be  very  important  for  our  new  long-term  effort  in  ionic  liquids.  The  derivation  of  an 
expression  for  the  analytic  gradient  for  the  EFP-ofe  initio  interaction  term  is  in  progress. 

A  very  important  new  development  (by  co-worker  Jan  Jensen)  is  a  new  method  for  using 
EFPs  across  covalent  bonds.  This  will  facilitate  the  representation  of  large  substituents,  as 
well  as  the  treatment  of  large  biological  molecules.  We  have  also  interfaced  the  EFP 
method  with  two  continuum  methods,  the  simple  Onsager  reaction  field  and  the  more 


sophisticated  polarizable  continuum  model  (PCM)10’11,  and  we  have  very  recently 
developed  a  density  functional  theory- based  EFP  method12. 

As  the  number  of  solvent  molecules  in  the  system  increases,  the  number  of  configurations 
to  be  considered  increases  rapidly,  and  traditional  small  molecule  geometry  optimization 
methods  are  not  effective.  We  are  therefore  developing  both  molecular  dynamics  and 
Monte  Carlo  simulation  codes  so  that  the  configurational  space  can  be  probed  more 
effectively,  not  only  for  minima,  but  for  transition  states  and  reaction  paths  as  well.  This 
is  also  an  important  first  step  in  the  development  of  methods  based  on  the  EFP  model  for 
predicting  bulk  properties  and  super-critical  behavior.  We  are  therefore  incorporating 
Ewald  summations  into  our  MD  method.  Collaborators  in  this  effort  are  Drs.  Paul  Day 
and  Ruth  Pachter  (AFRL)13.  Related  to  these  dynamical  methods  is  our  development  and 
implementation  of  a  method  (in  collaboration  with  Prof.  Michael  Collins,  Australian 
National  University)  for  converting  the  large  numbers  of  points  generated  in  ab  initio 
trajectories  into  a  global  potential  energy  surface,  using  a  modified  Shepard  interpolation 
approach. 

Major  James  Shoemaker’s  Ph.D  dissertation  in  Engineering  Physics  at  the  Air  Force 
Institute  of  Technology  focused  on  the  development  and  applications  of  an  embedded 
cluster  model  for  treating  surface  chemistry.  Papers  on  the  theoretical  method,  called 
SIMOMM,  has  been  published1'1’1".  Several  papers  applying  this  method  have  now  been 
published,  and  a  manuscript  that  describes  applications  to  silicon  carbide  surfaces  has 
been  published16.  Extensions  of  this  method  to  metal  oxide  catalysts  are  planned.  For 
metals,  however,  molecular  mechanics  is  not  likely  to  be  a  viable  approach  for  the  bulk, 
since  electrons  in  such  systems  are  too  delocalized.  To  overcome  this  problem,  we  have 
begun  to  explore  fast  multi  pole  methods  (FMM)  that  scale  linearly  and  are  also  highly 
parallel,  so  that  the  bulk  part  of  the  system  can  be  treated  by  quantum  mechanics17'18. 

Other  developments  .  We  have  also  been  exploring  alternatives  and  extensions  to  the 
popular  G2  and  G3  methods  developed  by  Pople  and  co-workers  for  the  accurate 
prediction  of  such  theormodynamic  properties  as  heats  of  formation,  ionization  potentials 
and  electron  affinities.  One  limitation  of  these  methods  is  that  they  are  applicable  only  to 
species  that  are  adequately  described  by  single  configuration  wavefunctions.  This 
eliminates  many  transition  metal  compounds,  di radicals  and  most  transition  states.  We 
have  therefore  developed,  in  collaboration  with  the  Radom  group,  multi- reference 
analogs  of  the  G2  and  G3  methods19.  The  multi -reference  methods  are  based  on 
CASSCF  wavefunctions,  instead  of  Hartree-Fock,  followed  by  multi-reference 
perturbation  theory  (instead  of  MJP2)  and  finally  multi-reference  Ci. 

Since  NMR  properties  are  a  very  important  experimental  diagnostic,  both  in  the  gas 
phase  and  in  solution,  we  are  in  the  process  of  developing  methods  and  codes  for 
predicting  NMR  properties.  Two  recent  accomplishments  are  the  derivation  for  the 
computation  of  NMR  coupling  constants  at  the  multi-reference  perturbation  level  of 
theory30  and  the  derivation,  implementation  and  first  successful  application  of  a  method 
to  calculate  NMR  chemical  shifts  in  the  presence  of  effective  fragment  potentials  . 
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B.  APPLICATIONS 


POSS  Compounds.  Polyhedral  oligomeric  silsesquioxanes  (POSS)  have  a  wide  variety  of 
important  applications,  including  lubricants  and  coatings.  They  are  also  resistant  to 
extreme  environments.  Consequently,  there  is  considerable  interest  in  these  species  in 
industry,  universities  and  Air  Force  Laboratories.  We  have  therefore  initiated  a  detailed 
study  of  POSS  compounds,  with  primary  emphasis  on  the  mechanisms  by  which  these 
species  form,  in  both  the  gas  phase,  and  in  solution  and  the  presence  of  catalysts.  In  the 
first  series  of  calculations,  we  examined  the  three  stepwise  hydrolysis  steps  of 
tricblorosilane,  followed  by  the  condensation  of  the  product  tri  hydroxy  silane.  All  four 
reactions  are  predicted  to  have  large  barrier  heights  in  the  gas  phase,  but  the  addition  of 
just  one  extra  water  molecule  is  sufficient  to  reduce  all  barriers  to  zero  or  nearly  zero, 
except  for  the  first  hydrolysis  step.  In  the  latter,  there  is  still  a  residual  barrier  of  almost 
10  kcal/mol.  So,  the  next  step  will  be  to  examine  the  effect  of  additional  water 
molecules,  especially  on  this  first  step. 

Calculations  have  also  been  completed  on  the  next  steps  in  the  mechanism,  in  which  the 
initial  condensation  products,  disiloxanes,  are  further  condensed  to  the  ring  compounds 
D3  and  D422.  These  results  are  consistent  with  those  reported  above,  in  that  initially  high 
barriers  are  reduced  to  nearly  zero  by  an  additional  water  molecule.  The  most  recent 
effort  has  been  the  study  of  the  addition  of  more  water  molecules,  which  lower  barriers 
even  further,  and  analysis  of  substituent  effects  on  the  reaction  mechanism  and 
energetics23. 

Although  our  main  focus  is  on  Si-based  POSS,  the  Ti  analogs  are  also  of  interest.  We 
have  therefore  initialed  a  series  of  calculations  on  the  Ti-POSS  compounds  to  study  their 
properties24.  These  mixed  Si-Ti  POSS  also  have  potential  as  catalysts,  and  a  paper  of  the 
efficacy  of  these  potential  catalysts  has  been  published25. 

Three  dimensional  cage  compounds,  including  zeolites,  have  been  of  interest  as  a 
possible  means  for  separating  small  gas  molecules.  In  collaboration  with  Dr.  Shawn 
Phillips  (AFRL-Ed wards)  we  have  therefore  completed  a  series  of  calculations  in  which 
the  potential  energy  surfaces  for  passing  N2  and  02  through  the  faces  of  T„,  n  =  8,  10, 

12  .  For  Tg,  both  the  energy  barriers  and  the  endothermicities  are  very  large  for  both  N2 
and  O2.  In  both  cases,  electron  correlation  is  essential  for  an  adequate  estimate  of  the 
energetics.  For  N2)  the  barrier  height  is  on  the  same  order  as  the  SiO  bond  strength, 
while  the  02  barrier  is  much  smaller.  Despite  these  observations,  the  SiO  bond  does  not 
break  upon  entry  of  either  molecule,  and  the  species  with  the  endohedral  gas  molecules 
are  minima  on  their  respective  potential  energy  surfaces.  As  the  size  of  the  cages 
increases,  the  barriers  and  endothermicities  decrease  as  one  would  expect. 

Fuels  with  High  N -Content.  There  has  been  considerable  interest  in  the  last  several  years 
in  the  potential  for  cubic  molecules  as  high  energy  fuels.  Cubane  itself  has  been 
considered,  and  so  far  rejected  because  of  the  complexity  of  its  synthesis.  The  octasila 
analog  of  cubane  has  been  prepared.  Unlike  the  carbon  system,  the  cubic  structure  is  the 
most  stable  SigHg  isomer.  This  makes  it  potentially  very  interesting,  since  it  is  still  rather 


high  in  energy.  However,  the  heavy  mass  of  Si  precludes  any  viability  of  such  species  as 
high-energy  fuels27. 

On  the  other  hand,  cubic  Ns  is  a  light,  very  high  energy  fuel  that  has  been  the  subject  of 
considerable  attention.  Extensive  calculations,  however,  illustrate  that  the  upper  limit  of 
the  barrier  separating  this  species  from  four  nitrogen  molecules  is  20  kcal/mol.  It  is, 
therefore,  not  a  viable  high  energy  fuel.  In  the  process  of  studying  the  cubic  Ns  potential 
energy  surface,  it  was  discovered  that  there  are  three  Ns  isomers  that  are  much  lower  in 
energy  than  the  cube,  yet  still  much  higher  in  energy  than  molecular  nitrogen.  The 
potential  energy  surfaces  for  these  isomers  have  now  been  studied  in  detail.  Two  of  these 
have  barriers  on  the  order  of  15-20  kcal/mol  (too  small),  but  the  third  has  a  barrier  of 
nearly  30  kcal/mol.  So,  this  isomer  may  be  a  viable  synthetic  target28. 

Related  to  the  foregoing,  much  attention  has  been  paid  recently  to  N4  as  a  potential  fuel. 
Many  references  are  made  to  a  barrier  on  the  order  of  50-65  kcal/mol,  even  though  it  has 
been  shown  by  Yarkony  that  non-adiabatic  interactions  reduce  this  barrier  to  less  than  40 
kcal/mol.  For  the  first  time,  we  have  mapped  out  the  potential  energy  surface  that  leads 
from  tetrahedral  N4,  using  second  order  perturbation  theory.  One  also  must  consider  the 
possibility  that  two  N4  molecules  could  interact  and  form  molecular  nitrogen.  This 
possibility  is  being  explored  in  this  laboratory. 

Christe  and  co-workers  at  AFRL  (Edwards)  have  isolated  the  first  new  all-nitrogen 
species  in  nearly  a  centuiy.  With  the  counter  ion,  this  species  is  [AsF6’][Ns+].  We  have 
now  studied  the  potential  energy  surface  of  this  species  as  it  forms  AsFs  +  FN5,  and 
subsequently  the  pathway  leading  from  FN5  to  FN3  +  N2.  A  paper  describing  this  work 
has  recently  been  published29.  Another  paper  on  the  possibility  of  making  crystals  of 
N5+/N5\  has  also  been  published30.  Both  of  these  latter  papers  have  been  done  in 
collaboration  with  colleagues  at  AFRL-Edwards. 

There  has  been  considerable  experimental  and  theoretical  interest  in  metal-doped  solid 
hydrogen.  We  have  an  ongoing  series  of  investigations  of  the  low-lying  electronic  states 
of  BH2  system,  to  determine  the  energetics  by  which  the  metal  is  held  as  a  weakly  bound 
species  in  the  matrix  and  the  ease  with  which  two  metal  atoms  find  each  other.  Similar 
studies  are  under  way  on  the  B2-H2  and  Al2-H2  systems. 

The  desired  reaction  product  for  burning  Al-doped  H2  is  AI2O3,  since  it  is  the 
thermodynamically  most  stable  aluminum  oxide.  However,  it  is  unknown  how  or  when 
this  species  forms  in  the  process.  We  have  therefore  initiated  an  extensive  study  of  the 
potential  energy  surface  that  leads  from  AI  and  02  to  various  oxides31  *32. 

II.  PROPOSED  RESEARCH 

A.  APPLICATIONS 

Ionic  Liquids.  We  will  initiate  a  broad  new  effort  in  the  field  of  energetic  ionic  liquids 
(EIL).  This  effort,  outlined  in  detail  in  the  following  paragraphs,  will  involve 
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collaborative  efforts  with  several  colleagues,  including  Drs.  Boatz,  Mills,  Drake,  and 
Hawkins  (AFRL-Edwards)  and  Professors  Thompson  (Oklahoma  State)  and  Voth  (Utah). 
Since  the  issues  involving  EIL  range  from  fundamental  questions  of  bonding,  kinetic 
stability  and  thermodynamics  to  condensed  phase  issues  such  as  phase  transitions, 
melting,  viscosity  and  surface  tension,  it  is  critical  that  we  interface  with  both  theoretical 
and  experimental  colleagues  in  order  to  make  optimal  progress. 

Requirements  for  a  viable  EIL  include: 

1 )  A  cation  that  has  a  large,  positive  heat  of  formation,  in  order  to  produce  as  large 
an  Isp  as  possible,  and  that  has  as  tow  a  melting  point  as  possible. 

2)  A  companion  anion  that  can  act  as  an  oxidizer. 

Therefore,  much  of  the  recent  focus  has  been  on  cyclic  nitrogen  containing  cations,  such 
as  triazoles  and  tetrazoles,  and  oxygen  containing  anions: 


1,2,3-triazoles 


1,2,4-triazoles 


Tetrazoles 
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Acyclic  cations  are  also  a  possibility,  but  the  strain  in  5-membered  rings  tends  to  increase 
the  heat  of  formation,  and  aromatic  character  can  stabilize  the  ring  sufficiently  to  make 
synthesis  feasible.  In  general,  asymmetry  in  the  ring  is  desirable,  since  this  seems  to 
lower  the  melting  point,  a  desirable  property.  Delocalization  of  the  ionic  charge  is  also 
desirable.  Target  anions  include  simple  ones  such  as  N03~,  ClOf,  and  N(NC>2)2 
(dinitramide),  as  well  as  more  complex  ones  like  NCIO3'2,  B(N03)4\  BfClO^f, 
A1(N03)4~.  Oxygen-free  anions  such  as  N3"  and  N5‘  are  also  of  interest. 

We  will  focus  on  both  sophisticated  electronic  structure  calculations  and  on  the  use  of  our 
effective  lfagment  potential  (EFP)  model  in  order  to  bridge  the  gap  between  the  gas  and 
condensed  phase.  The  former  will  be  in  collaboration  with  Dr.  Boatz,  while  the  latter 
will  rely  on  input  from  and  interactions  with  Professors  Thompson  and  Voth.  Specific 
efforts  are  summarized  in  the  following  projects.  Details  regarding  theoretical 
approaches  are  presented  in  Section  B. 

Properties  of  cation  monomers.  The  first  step  will  be  to  systematically  predict  the  G3 
proton  affinities  for  simple  amines,  R-NFF  +  Ef  ->  NH3R+.  This  simple  reaction  will 
provide  an  initial  indication  about  the  inherent  ability  of  a  particular  substituent  to 
increase  the  stability  of  the  amino  proton  and  therefore  resist  proton  transfer  to  the 
companion  anion.  Throughout  this  discussion,  the  R  groups  of  interest  will  include  N02, 
N3j  NH2,  CN,  CH3,  NF2,  CF3,  OH,  F,  C2H3,  OCH3.  The  next  step  will  be  to  optimize  the 
MP2  geometries  of  the  corresponding  substituted  triazoles  and  tetrazoles,  followed  by  the 
prediction  of  the  G3  heats  of  formation.  These  calculations  will  help  identify  those 
substituents  that  are  most  likely  to  increase  the  heat  of  formation  relative  to  the  parent 
compound.  For  those  substituents  that  are  the  most  promising,  the  corresponding  di¬ 
substitutions  will  be  investigated  as  well.  Although  the  main  interest  is  in  substitution  at 
N,  the  same  substituents  will  be  considered  at  C  as  well. 

Cation-anion  pair  interactions.  Once  the  inherent  heats  of  formation  as  a  function  of 
substitution  are  understood,  the  next  step  will  be  to  investigate  the  ion  pair  monomers. 
The  important  questions  regarding  these  simple  ion  pairs  are 

(1)  What  are  the  relative  energies  of  the  ion  pairs  relative  to  the  separated  ions? 

(2)  Is  there  significant  charge  transfer  upon  ion  pair  formation,  or  do  the  ions  retain  their 
distinct  ionic  character? 

(3)  For  those  species  in  which  proton  transfer  from  cation  to  anion  is  possible,  does  it 
occur  spontaneously?  If  not,  how  large  is  the  barrier  for  proton  transfer?  How  large  is 
the  barrier  for  transfer  of  the  R-(N)  group  from  cation  to  anion?  What  are  the  relative 
energies  of  ion  pairs  vs.  neutral  pairs? 

Monomer  ion  pairs  do  not  necessarily  reflect  the  behavior  of  ionic  liquids  in  the 
condensed  phase.  Therefore,  the  single  ion  pair  calculations  described  above  will  be 
foll  owed  by  dimer  calculations  to  determine  whether  the  presence  of  two  cation-anion 
pairs  serves  to  stabilize  the  species  by,  for  example,  inhibiting  proton  transfer  or  charge 


transfer.  Analogous  calculations  to  those  described  above  for  the  single  ion  pairs  will  be 
repeated  for  the  dimer  pairs.  Ultimately,  the  goal  is  to  build  a  unit  cell  for  each  of  the 
most  promising  species,  based  on  their  calculated  thermodynamic  properties.  The 
potential  energy  surfaces  of  these  species  will  provide  a  first  ( ab  initio)  step  in  the  study 
of  the  condensed  phases  of  these  systems.  The  latter,  discussed  in  more  detail  below,  will 
be  performed  in  collaboration  with  colleagues  (Thompson,  Voth)  who  are  expert  in  this 
field. 

All  of  the  geometry  optimizations  described  here  will  be  performed  using  second  order 
perturbation  theory  and  reliable  basis  sets.  For  the  smallest  species,  the  new  coupled 
cluster  codes  recently  implemented  in  GAMES S33  will  be  employed  to  test  the  MP2 
relative  energies.  Companion  calculations  will  be  performed  with  the  B3LYP  density 
functional  theory  approach,  in  order  to  assess  the  viability  of  this  approach.  The  B3LYP 
calculations  are  less  time-consuming  and  therefore  may  be  an  attractive  alternative  for 
larger  ion  pair  clusters. 

Decomposition  mechanisms.  All  of  the  foregoing  analyses  focus  on  structural  predictions 
and  thermodynamic  quantities  for  the  parent  compounds  and  ion  pairs.  The  next  step  will 
be  to  initiate  an  investigation  of  the  decompositions  of  the  triazole  and  tetrazole  rings. 
Since  it  is  not  known  whether  these  decompositions  occur  before  or  after  proton  (or 
substituent)  transfer,  both  possibilities  will  be  investigated.  Clearly,  it  is  not  possible  to 
know  beforehand  exactly  what  decomposition  mechanisms  will  most  viable,  however,  the 
following  general  approach  seems  sensible: 

(1)  Unimolecular  decompositions  of  the  cationic  and  neutral  ring  compounds  will  be 
studied.  Extrusion  of  a  very  stable  N2  moiety  seems  most  likely,  especially  for  the 
neutral,  unsubsituted  rings.  Substitution  at  various  positions  on  the  rings  is  likely  to  alter 
the  relative  stabilities  of  various  breakdown  products.  Barrier  heights,  minimum  energy 
paths  (MEPs),  and  ultimately  dynamic  reaction  paths  (DRPs)34  will  be  calculated  as  a 
function  of  substituent.  The  DRP  approach  (essentially  an  ab  initio  semi-classical 
trajectory  on-the-fly,  combined  with  an  RRKM  prediction  of  lifetimes,  was  particularly 
useful  in  our  analysis  of  the  FN5  decomposition  modes39.  A  DRP  can  be  performed  with 
any  level  of  theory  in  GAMES S,  and  the  availability  of  scalable  analytic  gradients  for 
MP2,  MCSCF,  and  DFT  makes  the  process  rather  efficient.  The  first  step  in  this  phase  of 
the  investigation  will  be  to  study  the  parent  ring  compounds.  This  will  provide  the 
baseline  against  which  the  effects  of  substituents  can  be  measured.  Then,  those 
substituted  compounds  that  are  most  promising,  based  on  the  earlier  thermodynamic 
studies  and  guidance  from  our  experimental  colleagues,  will  be  investigated. 

(2)  It  is  clearly  very  important  to  analyze  the  potential  energy  surfaces  for  the  oxidation 
of  the  cation  by  the  negatively  charged,  oxygen-containing  counter  ions.  Many  oxidation 
processes  are  possible  for  each  ring,  and  the  presence  of  substituents  may  very  well  alter 
the  mechanisms  dramatically.  As  for  the  unimolecular  decompositions,  the  first  step  will 
be  to  study  the  competing  oxidation  mechanisms  for  the  unsubstituted  compounds,  with 
an  emphasis  on  the  processes  that  are  likely  to  either  have  the  lowest  barriers  or  the  most 
stable  oxidation  products  (or  both).  Since  it  is  not  known  at  what  stage  the  oxidation 
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occurs,  both  the  ionic  and  neutral  species  wiil  be  studied,  using  the  same  computational 
methods  as  those  described  above.  By  this  point  in  the  project,  we  will  have  learned  a 
great  deal  about  the  effects  of  substituents  on  the  thermodynamic  and  unimolecular 
decomposition  processes.  Our  approach  to  the  oxidation  mechanisms  for  the  substituted 
compounds  will  be  guided  by  this  gained  understanding  and,  as  always,  by  our 
experimental  colleagues.  Since  there  is  currently  very  little  quantitative  experimental 
data  about  these  processes,  it  is  important  to  keep  in  mind  that  a  solid  qualitative 
overview,  especially  of  the  effects  of  substitution,  may  be  as  useful  as  detailed 
quantitative  comparisons. 

(3)  It  is  also  important  to  investigate  the  impact  of  heterogeneous  catalysis  on  the 
decomposition  mechanisms  and  energetics.  The  important  catalysts  are  Pt,  Pd,  Ni,  and  Ir. 
The  analysis  of  heterogeneous  catalysis  on  these  metals  will  be  accomplished  in 
increasingly  complex  stages.  Since  the  active  sites  for  catalysis  are  likely  to  be  relatively 
small  clusters  of  metals,  the  first  step  will  be  to  systematically  re-examine  the 
unimolecular  decomposition  and  oxidation  mechanisms  in  the  presence  of  small  metal 
clusters.  This  will  provide  us  with  an  initial  glimpse  into  the  manner  in  which  these 
mechanisms  are  altered  by  the  presence  of  the  catalyst.  We  are  currently  in  the  process 
of  developing  two-dimensional  periodic  boundary  conditions  for  the  study  of  surface 
science,  especially  heterogeneous  catalysis,  on  metal  surfaces  using  the  fast  multipole 
method  (FMM)  developed  already  for  both  the  Hartree-Fock  and  DFT  levels  of 
theory17’18,  and  extensions  of  this  approach  to  both  MCSCF  and  MP2  are  planned  (see 
Section  B  below).  Therefore,  in  the  longer  term  this  approach  will  facilitate  the  analysis 
of  the  heterogeneous  catalysis  of  the  decomposition  of  ionic  liquids  for  the  “infinite” 
systems. 

Extended  systems  and  connection  with  condensed  phase,  [  here  are  several  ways  in  which 
an  electronic  structure  theory  group  can  contribute  to  an  understanding  of  condensed 
phase  phenomena.  One  is  the  use  of  very  high  levels  of  theory  to  develop  accurate 
potential  energy  surfaces  that  can  then  be  used  by  condensed  phase  theorists  to  build 
reliable  potentials  with  which  to  perform  their  simulations.  The  proposed  calculations  that 
are  described  above  are  expected  to  serve  in  this  capacity,  in  addition  to  providing 
important  insights  in  their  own  right.  Second,  in  order  to  build  potentials  that  are  as 
reliable  as  possible,  it  is  also  important  to  develop  new  efficient,  highly  scalable  methods, 
so  that  the  most  sophisticated  methods  are  applicable  to  realistic  problems  on  tractable 
time  scales.  Several  such  proposed  developments  are  discussed  in  Section  B.  Third,  this 
group  has,  over  the  past  several  years,  been  developing  a  very  sophisticated  model 
potential,  the  effective  fragment  potential  (EFP)  that  is  firmly  based  in  quantum 
mechanics.  Described  in  detail  in  Section  B,  the  EFP  method,  when  completed,  wiil 
contain  all  of  the  fundamental  physics  involved  in  intermolecular,  and  inter- ionic, 
interactions.  Combined  with  periodic  boundary  conditions  and  either  FMM  or  Ewald 
sums  (both  of  which  will  be  coded),  the  EFP  method  will  allow  us  to  simulate  the  crystal 
phase  changes,  the  melting  process  and  subsequent  condensed  phase  behavior,  using  a 
combination  of  Monte  Carlo  and  molecular  dynamics  methods.  This  will  be  accomplished 
in  collaboration  with  Professor  Thompson.  The  EFP  capability  goes  well  beyond  that  of  a 
sophisticated  model  potential.  It  is  interfaced  with  the  fully  QM  capabilities  in  GAMESS, 
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so  that  it  can  also  function  as  a  QM/MM  method.  For  the  study  of  ionic  liquids,  this 
means  one  can  insert  a  QM  “solute”  such  as  a  single  ion  pair  (for  example)  into  the  “bulk” 
represented  by  EFPs  and  then  study  any  of  the  reaction  mechanisms  discussed  in  the 
preceding  paragraphs.  This  will  extend  our  understanding  of  these  systems  well  beyond 
the  small  clusters. 

Note  that  the  calculations  described  in  the  previous  paragraphs  will  be  very  demanding. 
Their  success  will  be  facilitated  by  the  continuing  development  of  the  parallel  code  in 
GAMESS  (See  Section  IIB)  and  continuing  access  to  the  DoD  Major  Shared  Resource 
Centers. 

Polynitrogen  species.  We  have  initiated  an  ongoing  collaboration  with  experimental 
colleagues  at  SRI  International  (the  late  Rob  Schmitt  and  Jeff  Bottaro)  to  investigate  the 
viability  of  a  number  of  polynitrogen  compounds  as  monopropellants.  The  genera! 
approach  is  to  first  estimate  the  heat  of  formation  of  each  species  using  variants  of  the 
G23=  and  G336  methods  developed  by  Pople  and  co-workers.  These  methods  make  use  of 
approximate  additivity  of  correlation  and  basis  set  corrections  to  obtain  thermodynamic 
properties  to  within  a  few  kcal/mol.  For  molecules  of  the  size  of  those  of  interest  here, 
these  are  very  demanding  calculations  and  require  the  use  of  parallel  computers  and  the 
parallel  electronic  structure  codes  developed  in  our  laboratory  (see  Section  B  below). 

Once  a  reliable  heat  of  formation  is  available,  colleagues  at  AFRL-Edwards  (Jerry  Boatz 
and  Jeff  Mills)  calculate  the  specific  impulse,  lsp.  If  the  Isp  seems  to  be  viable,  one  then 
must  consider  possible  decomposition  processes,  including  unimolecular  decomposition 
and  attack  by  obvious  environmental  species,  such  as  Nj,  O2,  and  water. 

The  initial  focus  in  this  phase  of  the  research  has  been  on  compound  1,  shown  below. 
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Using  a  combination  of  isodesmic  reactions  and  the  G2  model  on  the  IBM  SP2  at  the 
Maui  High  Performance  Computation  Center,  the  heat  of  formation  for  1  is  predicted  to 
be  456.8  kcal/mol.  This  translates  to  an  Isp  of  240  sec,  compared  with  230  sec  for 
hydrazine.  So.  this  appears  to  be  a  very  promising  fuel.  However,  it  is  always  important 
to  consider  the  stability  of  such  high  energy  species  to  various  reactions,  before  asserting 
their  viability  as  fuels.  One  possible  reaction  is  the  loss  of  molecular  nitrogen  from  the 
center  of  1  to  form  the  smaller  species  2  shown  below: 
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This  reaction  has  a  rather  large  60  kcal/mol  exothennicity.  The  barrier  height  for  this 
decomposition  must  be  determined,  since  given  the  complex  changes  involved,  could  be 
large  and  consequently  involve  a  large  barrier.  However,  1  could  also  break  an  N-N 
single  bond  at  either  or  both  ends  to  form  a  mono-  or  di -azide,  3: 

N3-C(N02)=N-N=N-N=C(N02)-M3  3 

This  may  still  be  an  interesting  species,  and  we  will  calculate  its  heat  of  formation,  but  it 
is  probably  considerably  less  energetic  than  1.  It  is  also  likely  that  loss  of  the  central  N2 
will  be  easier  in  3  than  in  f  and  equally  energetically  favorable.  Bottaro37  has  suggested 
that  4  (shown  below),  an  isomer  of  I,  may  be  more  viable  because  it  cannot  undergo  the 
analogous  NN  bond  cleavage  to  form  azides: 
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4 

It  is  likely  that  4  will  have  a  heat  of  formation  and  Isp  that  are  comparable  to  those  for  1. 
The  following  calculations  will  be  performed  on  this  molecule:  (1)  To  be  consistent  with 
the  calculations  on  1,  geometry  optimization  will  be  performed  using  second  order 
perturbation  theory  (MP2)  and  the  cc-pVTZ  basis  set;  (2)  The  most  likely  breakdown 
products  will  be  identified  (see  below),  the  corresponding  geometries  optimized,  and  the 
thermodynamic  energy  differences  calculated;  (3)  For  all  of  the  exothermic 
decompositions,  the  corresponding  transition  states  will  be  identified,  the  minimum 
energy  reaction  paths  followed,  and  the  dynamics  examined;  (4)  Assuming  the 
thermodynamic  and  kinetic  calculations  are  optimistic,  the  heat  of  formation  will  be 
calculated  in  the  manner  described  above  and  the  lsp  calculated.  The  simplest,  potentially 
exothermic  decomposition  one  can  imagine  for  4  is  the  extrusion  of  N2  from  one  or  both 
of  the  5-membered  rings.  Elimination  of  2  N2  would  yield  the  unusual  structure  4a 

02N-C-(N)6-C-N02  (4a) 

One  would  expect  4a  to  be  unstable  to  loss  of  two  additional  N2  molecules,  giving  two 
nitronitriles,  02N-C=N.  An  alternative  decomposition  route  for  4  would  be  to  eliminate 
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successive  O2N-ON  molecules,  ultimately  yielding  N&  which  will  then  rapidly  dissociate 
to  3  N2. 

Although  neither  1  nor  4  have  yet  been  synthesized,  the  related  compound,  5  below,  has 
been  made,  and  its  crystal  structure  has  been  determined  by  Gilardi  and  Karle38. 


5 

Using  approximate  additivity  relationships,  Bottaro  and  co-workers  have  estimated  the 
heat  of  formation  for  5  to  be  approximately  150  kcal/mol.  Since  this  is  a  known 
compound,  and  since  it  is  closely  related  to  compounds  I  and  4,  it  is  important  to 
determine  a  more  accurate  heat  of  formation  and  lsp  using  the  methods  described  above. 
We  will  also  examine  the  effects  of  replacing  the  hydrogens  with  amino  groups,  since  it 
is  estimated  that  each  such  group  will  add  about  47  kcal/mol  to  the  heat  of  formation37. 
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There  are  several  additional  N-0  compounds  that  the  SR!  group  have  targeted  for 
synthesis.  In  collaboration  with  Dr.  Jerry  Boatz  at  AFRL-Edwards,  we  will 
systematically  analyze  the  molecular  and  electronic  structures,  mechanisms  of  formation, 
and  heats  of  formation  for  these  compounds.  One  of  the  most  fascinating  of  these 
compounds  is  the  proposed  dimer  of  nitrosyl  azide,  N^NO.  The  structure  frequently 
written37  for  the  dimer  is  structure  (a)  below: 


However,  it  is  difficult  to  assign  electrons  using  this  arrangement  in  any  sensible  manner 
using  standard  Lewis  structures.  An  alternative  arrangement  is  the  eight-membered  ring, 
structure  (b),  in  which  there  are  alternating  single  and  double  NN  bonds  and  forma! 
single  bonds  to  the  oxygens.  This  structure  does  have  a  sensible  Lewis  structure  and  is 
also  rather  interesting  since  it  contains  an  eight-membered  nitrogen  ring.  This  ring  is 
isoelectronic  with  the  well-known  hydrocarbon  cyclooctatetraene,  so  it  is  likely  to  be 
nonplanar.  It  will  be  very  interesting  to  determine  the  actual  structure  of  this  dimer,  its 
stability  relative  to  the  monomers,  and  possible  pathways  (in  ground  or  excited  electronic 
states)  to  dimerization  and  decomposition.  As  always,  since  N2  is  so  stable,  elimination 
of  one  or  more  N2  molecules  seems  to  be  the  most  likely  decomposition  alternative  to 
simple  reversion  to  monomers.  Initial  calculations  on  these  species  will  be  performed 
using  MP2/cc-pVTZ,  however,  close  attention  will  be  paid  to  the  MP2  natural  orbital 
occupation  numbers,  in  order  to  monitor  possible  multi-reference  character.  The  unusual 
electronic  structure  of  (a)  suggests  significant  configurational  mixing  in  this  species,  not 
at  all  unusual  in  very  high  energy,  metastable  compounds.  These  compounds  will  also  be 
excellent  test  cases  for  the  new  completely  renormalized  (CR)-CCSD(T)  method  in 
GAM  ESS33.  The  CR-CCSD(T)  method  can  often  account  for  considerable  diradical 
character  when  normal  CCSD(T),  often  thought  of  as  state-of-the-art,  cannot. 

The  dimer  of  nitrosyl  azide  is  jut  one  example  of  a  series  of  N-oxides  which  the  SRI  team 
is  attempting  to  synthesize.  The  philosophy  here  is  that  the  presence  of  the  oxygen  atoms 
may  minimize  nitrogen  lone  pair  repulsions.  Additional  examples  are  the  two  cyclic 
anions  shown  below.  The  eight-membered  ring  structure  is,  of  course,  closely  related  to 
the  N3NO  dimer,  and  one  may  consider  both  the  neutral  and  dianionic  species. 
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Two  other  NO  anionic  species  whose  syntheses  are  actively  being  pursued,  and  which  we 
will  characterize  using  the  methods  discussed  above,  are  diaza-  and  triazanitrate: 


o  WN02) 

ii  I 

o2n-n-n-n-no2  o2n-n-n-n-no2 

It  is  important  to  emphasize  that  we  will  work  closely  with  the  SRI  group,  so  that  new 
calculations  can  be  initiated  as  new  promising  targets  are  identified. 

As  noted  in  Section  IA,  we  have  worked  with  the  AFRL-Edwards  group  to  characterize 
the  potential  energy  surface  for  [AsF6]'[N;>]+.  Christe  et  al.  have  prepared  [Sb2Fn]'[N5]+, 
and  the  new  counterion  apparently  considerably  increases  the  stability  of  the  [Ns]" 
species.  An  important  point  therefore  must  be  the  relative  net  energy  requirements  for 
the  F  transfer  to  [Ns]+  from  the  two  counterions.  We  will  use  second  order  perturbation 
theory  geometry  optimizations  to  determine  the  appropriate  energy  barriers  and  overall 
energy  changes  associated  with  the  transfer.  Ultimately,  one  clearly  must  understand  the 
solid  state  properties  of  these  compounds  in  order  to  fully  understand  their  stabilities. 

We  therefore  plan  to  analyze  the  potential  energy  surfaces  for  small  clusters  of  X'[N5]+, 
X  =  [AsF6]\  [Sb2Fn]\  in  order  to  enable  studies  by  colleagues  with  experience  in 
condensed  matter  theory  (e.g.,  Thompson,  Oklahoma  State;  Voth,  Utah). 

B.  THEORY  AND  CODE  DEVELOPMENTS 

Nearly  all  of  the  QM  calculations  described  above  will  require  highly  correlated  wave 
functions  and  reliable  atomic  basis  sets,  as  well  as  tools  to  study  environmental  effects. 
Because  of  the  complexity  and  size  of  most  of  the  species  of  interest,  it  will  be  important 
as  well  to  maximize  the  efficiency  of  the  calculations  with  highly  scalable  codes.  Some 
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of  the  methodologies  required,  most  notably  reliable  basis  sets,  are  being  developed  by 
others  and  will  simply  be  used  by  us.  For  example,  in  most  cases  the  basis  sets  of  choice 
will  be  those  collectively  referred  to  as  correlation-consistent  (cc)  basis  sets  developed  by 
Dunning  and  co-workers39.  Single  reference  coupled  cluster  methods  for  dosed  shells 
have  recently  been  added  to  GAMESS  (General  Atomic  and  Molecular  Electronic 
Structure  System)  in  collaboration  with  the  Piecuch  group  (Michigan  State)33.  These  new 
codes  also  include  the  completely  renormalized  (CR)-CCSD(T)  method  that  has  been 
shown  to  treat  bond-breaking  and  diradical  character  much  more  reliably  than  the  more 
traditional  methods.  Additionally,  in  collaboration  with  the  Piecuch  group,  we  are  about 
to  implement  the  equations-of-motion  (EOM)  coupled  cluster  methods  for  excited  states, 
including  the  EOM-CCSD(T)  method  that  is  not  readily  available  in  other  codes.  The 
next  steps  in  these  developments  by  the  Piecuch  group  will  be  to  implement  the 
analogous  open  shell  codes  and  to  derive  and  code  the  analytic  gradients.  Analytic 
gradients  for  a  given  CC  method  is  the  essential  step  in  developing  the  density  for  that 
method  and  therefore  all  of  the  interesting  response  properties  that  are  associated  with  the 
density.  So,  these  new  codes  that  are  being  developed  by  the  Piecuch  group  and 
implemented  into  GAMESS  by  the  ISU  group  will  add  significant  new  functionality  to 
GAMESS.  In  the  following  paragraphs  we  focus  on  those  new  theory  and  code 
developments  that  will  be  carried  out  at  Iowa  State  and  that  directly  impact  the  HEDM 
program. 

The  Effective  Fragment  Potential:  Solvent  Effects  and  Liquid  Properties.  Solvent  and 
other  environmental  effects  are  important  for  many  of  the  applications  discussed  above, 
in  particular  the  study  of  ionic  liquids.  In  addition,  there  are  many  problems  of  current 
interest  to  the  Air  Force  that  require  an  understanding  of  the  liquid-surface  interface.  In 
collaboration  with  Prof.  Jensen  (U.  Iowa),  we  will  continue  the  development  and 
implementation  of  new  techniques  to  treat  solvent  effects  and  liquids,  based  on  the 
effective  fragment  potential  (EFP)  method9,40,41.  Indeed,  the  ultimate  goal  of  the 
development  of  the  suite  of  EFP  methods  is  to  construct  model  potentials  from  first 
principles  that  are  able  to  address  the  broad  range  of  environmental  effects. 

Our  general  approach  to  solvent  effects  consists  of  three  layers.  The  innermost  layer 
consists  of  a  fully  quantum  mechanical  (QM:  ab  initio  or  DFT)  treatment  of  the  “solute” 
plus  n  solvent  molecules  (e.g.,  waters),  where  n  will  typically  be  a  small  number,  perhaps 
determined  by  the  role  particular  solvent  molecules  might  play  in  making  or  breaking 
new  bonds.  The  solute  may  be  a  single  molecule  of  interest  or  a  reacting  system.  The 
latter  is  important,  since  we  plan  to  study  environmental  (e.g.,  condensed  phase)  effects 
on  reaction  mechanisms  and  dynamics.  The  second  layer  consists  of  EFP  molecules:  a 
collection  of  m  explicit  solvent  molecules  in  the  immediate  vicinity  of  the  solute 
(compound  or  reacting  system).  The  internal  fragment  geometries  are  currently  held 
fixed,  but  the  positions  of  these  solvent  molecules  may  be  optimized.  One  may  think  of 
m  +  n  as  the  inner  solvation  layer.  The  third,  outermost  layer  is  some  representation  of 
the  bulk  (continuum)  solvent. 

The  interactions  between  solute  (QM)  and  solvent  (EFP)  molecules  are  treated  by  adding 
one-electron  terms  to  the  QM  Hamiltonian  to  represent  these  interactions9,40,41.  For  EFP- 
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EFP  interactions,  an  analogous  set  of  interaction  energies  are  derived.  For  each  QM-EFP 
and  EFP-EFP  interaction  term,  energy  gradients  are  derived  and  coded  to  enable 
geometry  optimizations,  reaction  path  following  and  dynamics  calculations  in  the 
presence  of  solvent  or  indeed  for  the  solvent  (liquid)  alone.  In  general,  such  a  potential 
may  be  expressed  as 


Eint  Ecoul  "t"  Elnd  F  Egr  "I"  Ect  +  Edisp  F  E]10, 

where  these  terms  represent,  respectively,  Coulomb,  induction  (polarization),  exchange 
repulsion,  charge  transfer,  dispersion  and  higher  order  interactions.  The  Hartree-Fock, 
GVB  or  MCSCF  levels  of  theory  contain  the  first  four  terms  that  are  determined  as 
follows9: 

-  Eeoui  is  treated  by  a  distributed  multipolar  expansion42  through  octupoles,  with 
expansion  points  at  the  atom  centers  and  bond  midpoints.  The  entire  expansion  is 
damped  by  a  screening  term  that  accounts  for  the  effect  of  overlapping  charge  densities 
as  the  distance  between  molecules  decreases. 

-  Eind  is  treated  in  a  self-consistent  manner  using  tensor  sums  of  localized  orbital 
polarizabilities  located  at  the  centroid  of  each  bond  and  lone  pair  of  the  fragment.  Both 
Ecoui  and  Eind  are  determined  from  QM  monomer  calculations. 

In  the  EFP1  method,  once  E^ui  +  Eind  have  been  determined  for  a  given  monomer  (e.g., 
H2O),  Eer  +  E(t  are  obtained  by  a  fitting  procedure.  One  determines  the  SCF  dimer 
interaction  energy  for  a  large  number  of  geometries.  For  water,  ~200  dimer  points  were 
obtained  by  choosing  several  water- water  orientations  and  then  several  0-0  distances  for 
each  orientation.  At  each  dimer  geometry  Ecou|  +  E,nd  is  subtracted  from  the  total  SCF’ 
interaction  energy,  and  the  remainder  (Eer  +  Ect)  is  fitted  to  gaussian  (EFP-EFP)  or 
exponential  (QM-EFP)  functions  that  are  centered  at  the  atom  centers  and  the  center  of 
mass.  Note  that  this  formulation  of  the  EFP  method  permits  systematic  improvements  in 
the  accuracy,  by  the  addition  of  higher  order,  but  still  easily  calculated  one-electron, 
terms  to  the  QM  Hamiltonian.  The  EFP  method  for  water  reproduces  ELF  results  with 
high  accuracy  for  geometries,  vibrational  frequencies,  energy  differences,  barrier  heights, 
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and  even  reaction  paths  .  In  order  to  incorporate  some  correlation  into  the  method,  a 
DFT-based  EFP1  method  has  recently  been  implemented12  using  the  B3LYP  functional. 
The  fitted  term  now  presumably  contains  some  short  range  correlation  in  addition  to  Eer  + 
Ect,  and  of  course  Ecoui  +  Emd  are  now  obtained  with  DFT. 

The  next  step  in  the  development  of  EFP  1  for  water  will  be  to  extend  the  method  to  the 
MP2  level  of  theory,  in  order  to  fully  incorporate  dispersion  effects.  EC0LI!  +  Ellld  will  be 
determined  in  a  manner  similar  to  that  described  above  for  SCF  and  DFT,  except  making 
use  of  the  MP2  multipoles  and  polarizabilities.  The  fitting  procedure  will  be  split  into  two 
pieces.  The  first,  again  representing  Ecoui  +  Emd  will  be  done  as  before  using  the  HF  fit. 
Then  the  sum  Ecoui  +  Emd  +  Eer  +  Ect  will  be  subtracted  from  the  each  of -200  MP2 
(H2CO2  interaction  energies.  This  remainder,  largely  Edisp,  will  then  be  fit  separately  to 
obtain  the  dispersion  function,  multiplied  by  a  damping  term  that  ensures  Edisp  0  at 
short  distances  (e.g.,  l-e'uR).  The  simplest  approach  is  to  assume  that  one  only  needs  the 
Q  term  in  the  dispersion  expansion42,  however,  we  will  explore  the  possibility  that 


higher  order  terms  (e.g.,  Cg)  might  improve  the  fit  to  MP2  interaction  energies.  The 
introduction  of  Edisp  and  the  corresponding  gradient  expression  will  also  require  the 
derivation  and  implementation  of  new  integrals  in  the  EFP  code,  since  the  Rys 
polynomials  will  not  be  applicable.  However,  these  will  still  be  simple  one-electron 
integrals. 

While  the  determination  of  the  E^ui  and  E^  are  easily  automated  and  need  only  be  done 
once  for  a  given  solvent,  the  fitting  process  for  Eer  +  ECI  and  E<jisp  is  a  significant 
bottleneck,  since  it  must  be  repeated  for  each  solvent  of  interest.  For  more  complex 
solvents  the  number  of  necessary  dimer  points  may  be  very  large.  Therefore,  in  addition 
to  the  EFP1  method  for  H20,  we  have  begun  to  explore  a  more  general,  “first  principles” 
EFP2  method  that  does  not  rely  on  fits  to  QM  data.  This  has  already  been  accomplished 
for  the  exchange  repul  sion  between  two  fragments*1,  by  making  use  of  an  expansion  in 
the  overlap  matrix,  terminated  at  second  order.  By  representing  the  wavefunctions  on  the 
interacting  molecules  by  localized  molecular  orbitals,  this  simple  method,  that  requires 
only  the  evaluation  of  one-electron  integrals,  provides  a  veiy  accurate  representation  of 
the  exact  intermolecular  exchange  energy.  This  method  has  now  been  incorporated  into 
the  EFP  code,  complete  with  analytic  gradients,  and  it  has  been  successfully  tested  for 
many  solvents9.  The  extension  of  this  method  to  the  QM-EFP  exchange  repulsion  has 
also  been  completed  and  incorporated  into  GAMESS. 

So,  the  current  EFP2  method  includes  Eer  +  Ect  +Eer  with  no  fitted  parameters.  The  next 
steps  will  be  to  first  derive  and  code  the  analytic  gradient  for  the  QM-EFP  Eer  interaction. 
Second,  we  will  develop  an  independent  term  for  Edisp-  Dispersion  will  be  introduced  into 
EFP2  initially  by  following  the  lead  of  Amos  et  al45  ,  by  developing  and  implementing  a 
distributed  LMO  model, 
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where,  for  example,  ak(iv)  is  the  frequency-dependent  polarizability  of  LMO  k  at 
imaginary  frequency  v  and  7^  is  the  field  gradient  (which  scales  as  r 3).  The  integral  will 
be  computed  by  numerical  quadrature  using  ~20  points  as  suggested  by  Stone.  Thus, 
the  frequency-dependent  polarizability  tensors  will  be  pre-computed  for  each  LMO  at 
~20  frequencies  using  time-dependent  (TD)  DFT48’47  and  stored  as  part  of  the  EFP. 
a\iv)  tends  to  zero  with  increasing  frequency,  so  the  proper  frequency  range  can  easily 
be  determined  by  systematically  decreasing  the  frequency.  TD  DFT  is  expected47  to 
predict  a(iv)  to  an  accuracy  of  10-20%.  Since  the  a(iv)  need  only  be  calculated  once 
for  a  given  molecule,  this  approach  can  be  tested  for  small  molecules  (e.g.,  H20,  NH3) 
using  our  new  fully  parallel  Full  Cl  code48  and  the  sum  over  states  (SOS)  method.  This 
will  provide  exact  results  for  a  given  atomic  basis  and  will  therefore  calibrate  the  TD 
DFT  method  that  must  be  used  for  larger  molecules  of  interest.  Since  charge  penetration 
effects  can  significantly  damp  the  dispersion  energy42  at  short  range,  we  will  initially 
explore  the  damping  approach  introduced  by  Thole,49  in  which  the  field  gradient  is 
damped  by  modifying  its  functional  form  to  reflect  an  approximate  charge  distribution 
shape.  While  Thole’s  approach  introduces  empirical  parameters,  we  will  explore  the  use 
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of  the  SGO  (spherical  Gaussian  orbital)  approximation9,44  to  introduce  overlap-dependent 
parameters  that  can  be  evaluated  directly. 


The  QM-EFP  dispersion  energy  analogous  to  Eq  (I )  can  be  derived  by  a  multipole 
expansion  of  EFP  B  only, 
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The  integral  is  over  orbitals  on  QM  atom  A  and  S;  is  the  energy  of  orbital  i.  Thus,  the 
dispersion  energy  is  computed  using  converged  SCF  orbitals  in  analogy  with  the  MP2 
energy.  Since  both  equations  only  involve  one-electron  integrals  they  can  be  computed  at 
trivial  expense  compared  to  the  MP2  energy  of  molecule  A.  The  damping  contribution 
will  be  incorporated  by  modifying  the  electric  field  operator  as  described  above.  The 
energy  and  gradient  expressions  will  be  derived,  coded  and  tested  as  part  of  the  proposed 
research. 
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Third,  Jensen50  has  derived  an  expression  for  Ect,  using  the  same  SGO  approach  that  has 
been  so  successful  for  Eer.  The  gradient  for  this  term  will  be  derived,  and  the  energy  + 
gradient  tested  for  accuracy  and  efficiency.  A  proper  treatment  of  Ect  will  be  especially 
important  for  ionic  interactions.  So,  this  will  be  a  key  development  for  our  ability  to  treat 
ionic  liquids. 


The  EFP  method  has  already  been  interfaced  with  MCSCF  wavefunctions  in  GAMES S, 
so  qualitative  studies  of  excited  state  solvent  effects  are  already  possible51.  Nonetheless, 
an  interface  with  the  Cl  code  is  desirable  for  more  accurate  excited  state  solvation 
studies.  Webb52  has  very  recently  implemented  a  distributed  parallel  Cl  singles  method 
into  GAMESS,  including  analytic  gradients.  The  gradients  provide  us  with  the 
Hellmann-Feynman  term  for  the  Cl  first  order  density  matrix,  so  we  are  now  in  a  position 
to  develop  a  Cl- EFP  interface  for  treating  excited  state  solvent  effects,  including 
geometries,  reaction  paths,  and  dynamics.  These  capabilities  will  be  implemented  during 
the  next  grant  period.  A  longer  term  goal  is  to  derive  and  implement  gradients  for  multi¬ 
reference  Cl  and  EOM-CC  methods,  so  that  their  first  order  density  matrices  will  be 
available  to  interface  with  EFP2. 


In  order  to  treat  extended  systems,  such  as  those  to  be  encountered  in  the  ionic  liquids 
project,  it  will  be  necessary  to  develop  both  Monte  Carlo  and  molecular  dynamics 
methods  with  periodic  boundary  conditions  and  long-range  cutoffs.  A  Monte  Carlo  (with 
basin  hopping)  simulated  annealing  code  has  been  implemented  and  tested  on  water 
clusters13  and  simple  organic  acids53.  For  the  study  of  ionic  liquids,  we  will  need  to 
extend  this  code  to  include  periodic  boundary  conditions,  and  either  Ewald  sums  or  FMM 
for  long-range  cutoffs.  Since  we  have  already  developed  FMM  methods17,18,  this  will  be 
the  first  approach,  however,  Ewald  sums  will  also  be  coded.  This  new  MC  code  will  then 
facilitate  the  study  of  phase  changes  and  melting  of  ionic  liquids.  We  have  also 
implemented  an  EFP1/MD  code,  using  periodic  boundary  conditions  (PBC),  the 
minimum  image  convention  (MIC)  and  quaternions  to  deal  with  the  internally  frozen 
geometries.  A  fully  QM  MD  code  is  also  in  progress.  As  part  of  the  proposed  research, 
we  will  extend  this  code  to  the  new  EFP2  method,  complete  with  dispersion  and  charge 


20 


transfer  terms,  and  then  implement  Ewald  sums  to  deal  with  long-range  interactions.  The 
combination  of  MC  and  MD  methods  will  allow  us  to  simulate  bulk,  condensed  phase 
behavior.  Once  the  new  QM/EFP  interaction  terms  (dispersion,  charge  transfer)  have 
been  derived  and  coded,  the  Monte  Carlo  code  will  be  extended  to  QM/EFP2 
calculations.  Likewise,  once  analytic  gradients  are  derived  and  coded  for  each  of  the 
QM/EFP2  interaction  terms,  the  MD  code  will  be  adapted  to  QM/EFP2  calculations. 

This  combination  of  techniques  will  facilitate  the  study  of  the  decomposition  mechanisms 
of  ionic  liquids  in  the  presence  of  the  extended  system. 

Scalable  Computing  and  Linear  Scaling  Methods.  The  computer  time  requirements  for 
ab  initio  calculations  increase  dramatically  with  the  size  of  the  system  to  be  treated,  even 
at  the  SCF  level.  Even  model  potentials,  which  require  orders  of  magnitude  less 
computation  time  than  ab  initio  computations  for  an  energy  +  gradient  run,  can  be 
extremely  time-consuming  if  one  requires  tens  of  thousands  of  points  in  a  molecular 
dynamics  (MD)  or  Monte  Carlo  simulation.  Two  very  successful  approaches  for  solving 
this  problem  in  our  group  have  been  the  development  of  new  methods  for  scalable, 
distributed  parallel  computing  and  new  approaches  for  linear  scaling  algorithms.  Much 
of  the  parallel  developments  in  GAMESS34  over  the  past  several  years  has  been 
supported  by  a  separate  DOD  CHSSI  software  development  grant  that  will  almost 
certainly  not  be  renewed  beyond  the  termination  date  of  12/04.  That  grant  has  enabled  us 
to  develop  the  distributed  data  interface  (DDI)33’56  and  thereby  implement  algorithms  for 
MP2  and  UMP2  energies  +  gradients,  and  ZAPT2  and  MRMP2  energies,  that  scale 
essentially  perfectly  up  to  5 12  processors  on  systems  like  the  T3E  and  to  32-64 
processors  on  UNIX  or  Linux  clusters.  A  new  DDl-based  scalable  CASSCF  analytic 
Hessian  code  has  just  been  completed37  and  appears  to  scale  very  well.  Analogous  codes 
are  under  development  for  ZAPT2  gradients  and  CASSCF  energies  +  gradients. 

The  main  focus  regarding  parallel  code  development  for  the  research  proposed  here  will 
be  on  scalable  QM/EFP  methods  discussed  above,  and  the  linear  scaling  methods 
discussed  below.  The  over-riding  motivation  is  to  enable  us  to  use  larger  QM  regions  in 
the  QM/EFP  calculations,  and  to  perform  simulations  using  larger  EFP  regions.  The  one 
addition  is  the  need  to  develop  scalable  gradients  for  state-averaged  CASSCF 
calculations,  because  this  method  is  essential  for  the  calculation  of  non-adiabatic  matrix 
elements,  such  as  spin-orbit  coupling  and  vibronic  (derivative)  coupling.  The  latter  could 
well  be  important  for  the  study  of  some  of  the  potential  energy  surfaces  of  poly  nitrogen 
species. 

As  part  of  the  CHSSI  program,  we  have  developed  codes  for  EFP1/HF  that  scale  very 
well  to  32  processors,  as  long  as  one  sensibly  grows  the  size  of  the  system  with  the 
number  of  processors.  This  effort  benefited  greatly  from  the  efforts  of  Heather  Netzloff, 
a  student  who  was  supported  entirely  by  her  own  fellowship  funds  throughout  her 
graduate  career.  The  analogous  code  for  EFPI/DFT  should  now  be  straightforward,  and 
this  will  be  our  first  goal.  The  parallel  code  for  EFP1/MP2  should  be  similar,  as  the  only 
completely  new  terms  will  be  those  for  the  dispersion  energies  and  gradients.  Once 
parallel  codes  for  all  of  the  EFP1  variants  have  been  completed,  a  parallel  code  for  the 
EFP2  method  will  be  developed.  Since  the  Coulomb  and  polarizability  terms  will  be 
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analogous  to  those  in  EFP1,  the  focus  will  naturally  be  on  the  exchange  repulsion, 
dispersion  and  charge  transfer  energies  and  corresponding  gradients. 

The  discussion  in  the  previous  paragraph  focused  entirely  on  developing  parallel  EFP 
codes,  but  developing  the  analogous  codes  for  the  EFP-QM  interaction  terms  is  at  least  as 
important.  One  would  like,  for  example,  to  have  the  capability  to  perform  Monte  Carlo 
or  molecular  dynamics  simulations  on  EFP  arrangements  in  the  presence  of  a  QM 
“solute”.  In  order  to  make  such  a  calculation  feasible  for  a  reasonable  number  of 
structures,  one  needs  all  EFP-EFP  and  EFP-QM  terms  to  be  calculated  as  efficiently  as 
possible.  Given  our  success  in  developing  such  codes  for  QM  and  EFP  separately,  there 
is  little  doubt  that  this  will  be  a  successful  effort.  Several  of  the  QM  methods  in 
GAMESS  (e.g.,  CCSD(T),  MRMP2)  do  not  yet  have  analytic  gradients.  Fortunately,  we 
have  recently  developed  a  group-based  extension  of  DDI  (GDDI)38'*9  in  which  numerical 
derivatives  become  trivially  parallel  in  a  coarse-grained  sense  by  having  each  function 
evaluation  performed  on  a  different  node.  If  the  function  itself  exists  as  a  parallel  code  in 
GAMESS  (e.g.,  MRMP2),  and  if  the  nodes  are  SMP  nodes,  the  multi-level  parallelism 
can  be  exploited.  At  present,  GDDI  can  only  be  applied  to  fully  QM  calculations.  As 
part  of  the  proposed  research,  we  will  extend  this  very  useful  approach  to  EFP 
calculations. 

Building  on  beautiful  work  by  the  Head -Gordon60  and  other61  groups,  we  have  previously 
developed  an  optimized  linear  scaling  fast  multipole  (FMM)  method  for  the  Coulomb 
term  in  both  HF  and  DFT17,18.  The  novelty  of  this  approach  is  that  it  is  formulated  in 
such  a  manner  that  one  can  identify  a  small  number  of  parameters  that  can  be  easily 
optimized  in  an  iterative  manner  to  optimize  the  linearity  of  the  scaling.  The  desired 
accuracy  for  the  entire  calculation  is  an  input  parameter.  A  replicated  data  parallel  code 
was  developed  for  this  method  that  scales  moderately  (factor  of  six  on  eight  processors 
for  a  model  problem).  Several  new  features  will  be  developed  as  part  of  the  proposed 
research :  (1)  Analytic  derivatives  will  be  derived  and  coded,  to  enable  geometry 
predictions,  reaction  path  following  and  the  evaluation  of  vibrational  frequencies  via 
semi-numerical  Hessians;  (2)  Flexibility  will  be  introduced  into  the  method  so  that 
accuracy  in  some  regions  of  an  extended  system  (e.g.,  a  protein)  that  are  far  from  the 
chemical  “action”  are  allowed  to  be  lower  than  the  accuracy  required  where  the 
chemistry  is  taking  place;  (3)  The  method  will  be  extended  to  MCSCF  wavefunctions, 
since  many  problems  of  interest  to  us  are  inherently  multi-configurational,  and 
subsequently  to  second  order  perturbation  theory  methods;  (4)  Distributed  data  parallel 
codes  will  be  developed  for  all  levels  of  our  linear  scaling  ansatz.  All  of  these  techniques 
will  be  incorporated  into  the  2-D  PBC  code  we  are  already  developing  to  treat  the 
chemistry  that  occurs  on  metal  surfaces.  This  will  provide  the  essential  ingredients  for 
the  study  of  heterogeneous  catalysis  of  ionic  liquids  decomposition  on  metal  catalysts. 

Thermodynamic  properties.  The  availability  of  accurate  methods  for  predicting 
thrmodynamic  properties,  such  as  heats  of  formation,  is  important  in  order  to  accurately 
estimate  the  specific  impulse.  The  most  popular  such  methods  are  the  G23:i  and  G336 
methods  developed  by  Pople  and  co-workers.  One  limitation  of  these  methods  is  that 
they  are  applicable  only  to  species  that  are  adequately  described  by  single  configuration 
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wavefunctions.  This  is  a  serious  limitation  for  many  species  of  interest  to  the  HEDM 
program,  since  such  species  frequently  have  unusual  bonding  arrangements  and  are 
therefore  not  well  described  by  single  configuration  wavefunctions.  We  have  therefore 
developed,  in  collaboration  with  the  Radom  group,  multi-reference  analogs  of  the  G2  and 
G3  methods19.  In  these  methods,  the  single  configuration  Hartree-Fock  wavefunction  is 
replaced  by  a  full  valence  MCSCF  wavefunction,  single  reference  MP2  is  replaced  by 
multi-reference  second  order  perturbation  theory  (MRMP),  and  QC  ISD(T)  is  replaced  by 
multi-reference  Cl  (MRCI).  The  method  has  been  tested  on  the  G2-1  test  set,  with 
accuracy  similar  to  that  of  the  original  single  configuration-based  methods62. 
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